function err_mom = objective(x, Markuptarget)

optset('bisect', 'tol', 1e-18); 

% Calibrated Parameters

p.xi       = x(1);                        % Pareto tail productivity
p.sigma    = x(2);                        % demand elasticity: match level of markups 

% Assigned Parameters

p.phi      = 0.4;                         % doesn't matter since reset below to match material share  

p.epsi     = 0.162*p.sigma; 

p.beta     = 0.96;                        % period is 1 year
p.varphi   = 0.04;                        % exit rate 
p.delta    = 0.06;                        % depreciation rate

p.r        = 1/p.beta - 1; 
p.nu       = 1;                           % inverse labor supply elasticity
p.alpha    = 1/3;                         % capital share in value added
p.theta    = 0.5;                         % elasticity of substitution between VA and intermediates

p.xie      = 0;                           % entry subsidy
p.xis      = 0;                           % sales subsidy
p.taus     = 0;                           % sales tax (to isolate role of removing misallocation) -- shows up in planner's problem

% Initial Aggregates (Normalization)

mshare     = 0.45;                      % share of intermediate inputs in total sales

Y          = 1; 
N          = 1; 
R          = p.r + p.delta; 

printr     = 0; 

if printr == 1
    
    options.Display = 'iter';

else
    
    options.Display = 'none';

end


% Parameters governing accuracy of approximation

% Quadrature
 
nz         = 5000;                         % Gauss-Legendre quadrature for exponential r.v.

zmax       = -1/p.xi*log(1e-22);

[z, w]     = qnwlege(nz, 0, zmax); 
w          = p.xi.*w.*exp(-p.xi.*z);
w          = w./sum(w); 

z          = exp(z); 

p.w        = w; 
p.z        = z; 


if printr
    
loglog(z, 1 - cumsum(w));
ylim([1e-8, 1])

fprintf('Mean z (True, Quadrature)     = %9.4f %9.4f\n',  [p.xi/(p.xi - 1),               w'*z]);
fprintf('Var z (True, Quadrature)      = %9.4f %9.4f\n',  [p.xi/(p.xi - 1)^2/(p.xi - 2),  w'*(z - w'*z).^2]);
fprintf('E z^4 (True, Quadrature)      = %9.4f %9.4f\n',  [(p.xi/(p.xi - 4))^(1/4),       (w'*z.^4)^(1/4)]);

end

%%%%% 1. Solve D and firm's choices for a given N %%%%%% 

%%%% Start solving

D              = fsolve('findequilibrium', 1, options, w, z, p, 'market', 'old');

[~, ~, ~, Omega, ~,~, ~, ~, ~, ~, ~, ~, Mu]  = findequilibrium(D, w, z, p, 'market', 'old');

p.phi          = 1 - mshare*Omega^(1 - p.theta)*Mu;                                       % choose to match materials share in sales = mshare 


[~, q, mu, Omega, Z, W, Lp, K, B, U, Uq, C, Mu]  = findequilibrium(D, w, z, p, 'market', 'old');

p.F            = Y/N/W/(1/p.beta - 1 + p.varphi)*(1 - 1/Mu); 

L              = Lp + p.F*p.delta*N;

p.psi          = W/(C*L^p.nu); 

P              = (N*w'*(Uq.*q))^(-1);       % Demand index so we can compute prices


% Individual Choices

py             = Uq*P;                      % producer price
y              = q; 
l              = py.*y/Y.*Mu./mu.*Lp; 
b              = py.*y/Y.*Mu./mu.*B; 
k              = py.*y/Y.*Mu./mu.*K; 
mc             = Omega./z; 

if printr
    
% Report Some Aggregate Statistics
    
fprintf('\n');
fprintf('\n')
fprintf('Entry cost, kappa                       = %9.3f \n',  p.F);
fprintf('Disutility from work, psi               = %9.3f \n',  p.psi);
fprintf('\n');
fprintf('\n');
fprintf('Measure firms, N                        = %9.3f \n',  N);
fprintf('Capital-Output rato, A/Y                = %9.3f \n',  K/Y);
fprintf('N-Output rato, N/Y                      = %9.3f \n',  N/Y);
fprintf('Wage rate                               = %9.3f \n',  W);
fprintf('Output, Y                               = %9.3f \n',  Y);
fprintf('Consumption, C                          = %9.3f \n',  C);
fprintf('Materials, B                            = %9.3f \n',  B);
fprintf('GDP, Y - B                              = %9.3f \n',  Y - B);
fprintf('Investment, X                           = %9.3f \n',  p.delta*K);
fprintf('Employment, L                           = %9.3f \n',  L);
fprintf('Employment (production), Lp             = %9.3f \n',  Lp);
fprintf('Labor share (production)                = %9.3f \n',  W*Lp/Y);
fprintf('Agg Labor share                         = %9.3f \n',  W*L/Y);
fprintf('Variable Cost share                     = %9.3f \n',  (W*Lp + B)/Y);
fprintf('\n');
fprintf('Investment to GDP                       = %9.4f \n',  p.delta*K/(Y - B));
fprintf('Materials Share in Sales                = %9.4f \n',  B/Y);
fprintf('Profits to GDP                          = %9.4f \n',  (Y - W*Lp - R*K - B)/(Y - B));

fprintf('\n');

weight = w.*l/sum(w.*l); 

fprintf('10 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 10, weight));
fprintf('25 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 25, weight));
fprintf('50 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 50, weight));
fprintf('75 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 75, weight));
fprintf('90 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 90, weight));
fprintf('95 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 95, weight));
fprintf('99 p.c Markup (cost-weighted)           = %9.3f \n',  wprctile(mu, 99, weight));


fprintf('\n');
fprintf('\n');

fprintf('Various Checks of the Code\n');
fprintf('\n');
fprintf('Error in Kimball aggregator             = %9.2e \n',  N*w'*U - 1);
fprintf('Zero Profits Final Goods Firm           = %9.2e \n',  N*w'*(py.*y) - Y);
fprintf('Error in Consumption-Leisure Choice     = %9.2e \n',  p.psi - W/(C*L^p.nu));
fprintf('Error in Capital Resource constraint    = %9.2e \n',  N*w'*k - K);
fprintf('Error in Labor Resource constraint      = %9.2e \n',  N*w'*l - Lp);
fprintf('Error in Materialsn Res constraint      = %9.2e \n',  N*w'*b - B);
fprintf('1 - Markup x Marginal Cost Aggreg       = %9.2e \n',  1 - Mu*Omega/Z);

fprintf('Error in Markups                        = %9.2e \n',  norm(py(q > 0)./mc(q > 0) - mu(q > 0)));

fprintf('\n');
fprintf('\n');


Dp       = fsolve('findequilibrium', D, options, w, z, p, 'planner', 'old');

[~, qp, ~, ~, Zp, ~,~, ~, ~, Up]  = findequilibrium(Dp, w, z, p, 'planner', 'old');

fprintf('\n');
fprintf('Error Kimball aggregator planner             = %9.2e \n',  N*w'*Up - 1);
fprintf('\n');
fprintf('Losses from Misallocation: Z, p.p.           = %9.2f \n',  log(Zp/Z)*100);

fprintf('\n');
fprintf('\n');

end

    
moment_data        = [Markuptarget; 0.45; 0.229; 0.573; 0.739; 0.162]; 



moment_model       = zeros(size(moment_data)); 

moment_model(1)    = Mu; 
moment_model(2)    = B/Y;


data               = sortrows([w, py.*y], 2); 

cumF               = cumsum(data(:,1)); 


top1               = cumF >= 0.99; 
top5               = cumF >= 0.95; 
top10              = cumF >= 0.90; 

moment_model(3)    = data(top1,  1)'* data(top1,  2)/(data(:,1)'*data(:,2));
moment_model(4)    = data(top5,  1)'* data(top5,  2)/(data(:,1)'*data(:,2));
moment_model(5)    = data(top10, 1)'* data(top10, 2)/(data(:,1)'*data(:,2));


keep               = py.*y > 1e-8; 

bhat               = lscov([ones(size(y(keep))), log(py(keep).*y(keep))], 1./mu(keep) + log(1 - 1./mu(keep)));

moment_model(6)    = bhat(2);


if printr
    

fprintf('\n');

fprintf('Moments used in Calibration: Left = Model, Right = Data  \n');   
fprintf('\n');
fprintf('\n');
fprintf('Aggregate Markup                           = %9.3f %9.3f\n',  [moment_model(1),    moment_data(1)]);
fprintf('Share Intermediates in Sales               = %9.3f %9.3f\n',  [moment_model(2),    moment_data(2)]);
fprintf('\n')
fprintf('Sales share largest 1   percent firms      = %9.3f %9.3f\n',  [moment_model(3),    moment_data(3)]);
fprintf('Sales share largest 5   percent firms      = %9.3f %9.3f\n',  [moment_model(4),    moment_data(4)]);
fprintf('Sales share largest 10  percent firms      = %9.3f %9.3f\n',  [moment_model(5),    moment_data(5)]);
fprintf('\n');
fprintf('Regression coefficient                     = %9.3f %9.3f\n',  [moment_model(6),    moment_data(6)]);
fprintf('\n');
end


weights        = zeros(numel(moment_model), 1); 

weights([1; 2; 4])  = 1; 

weights        = weights/sum(weights);


err_mom        = (moment_model - moment_data)./(1 + moment_data);
err_mom        = (weights'*err_mom.^2).^(1/2);

if printr
        
fprintf('\n');
fprintf('\n');   
fprintf('Error in Moments                           = %9.6f \n', err_mom);
fprintf('\n');
fprintf('\n');   

end

format short g
fprintf('\n');
disp([x(:)',  err_mom]);
